Reservoirs.f90 Source File

!! Simulate reservoirs within the distributed hydrological model
!|author:  <a href="mailto:giovanni.ravazzani@polimi.it">Giovanni Ravazzani</a>
! license: <a href="http://www.gnu.org/licenses/">GPL</a>
!    
!### History
!
! current version  2.0 - 4th March 2025  
!
! | version  |  date       |  comment |
! |----------|-------------|----------|
! | 1.0      | 05/Oct/2016 | Original code |
! | 1.1      | 08/Feb/2023 | OutReservoirsInit and OutReservoirs added |
! | 1.2      | 17/Mar/2023 | bypass channel simulation included |
! | 1.3      | 28/Nov/2023 | manage reservoir when level is high |
! | 1.4      | 12/Dec/2023 | eflow is a vector of 365 daily value |
! | 1.5      | 26/Mar/2024 | reservoir volume written to output |
! | 1.6      | 10/Apr/2024 | observed reservoir downstream discharge read from file |
! | 1.7      | 11/Apr/2024 | observed reservoir diverted discharge read from file |
! | 1.8      | 05/May/2024 | Qin and Qout channel written even in reservoir without diversion |
! | 1.9      | 04/Sep/2024 | modified code to save and read reservoir state variables |
! | 2.0      | 04/Mar/2025 | one outlet discharge rule can be assigned for every day of the year |
!
!### License  
! license: GNU GPL <http://www.gnu.org/licenses/>
!
!### Module Description 
! This module includes data and routines to manage 
!  reservoirs within the distributed hydrological model
! Two types of reservoirs are available, _on-stream_ and 
! _off-stream_ reservoirs. 
! _On-stream_ reservoirs are the ones that are located 
! across a river. An _off-stream_ reservoir is a reservoir 
! that is not located on a streambed, and is supplied by 
! an artificial canal or pipeline. In case of _off-stream_
! reservoir, the water stored can be released in a 
! downstream section of the same river from which water 
! was withdrawn, or in a different water course.
! A _diversion channel_ can be associated to a given
! reservoir. In this case, the _diversion channel_  
! diverts water from the reservoir and conveys it 
! to a downstream section of the same river (bypass) 
! or to a different river. The  _diversion channel_ 
! is an hydraulic structure typically used to
! diverts flow to hydroelectric power plant.
! 
MODULE Reservoirs      

! Modules used: 

USE DataTypeSizes, ONLY : &  
  ! Imported Parameters:
  float, short

USE Chronos, ONLY : &
  !Imported definitions:
  DateTime, &
  !Imported variables:
  timeString, &
  !Imported operators:
  ASSIGNMENT (=)
  
USE ObservationalNetworks, ONLY : &
  ! Imported definitions:
  ObservationalNetwork, &
  !Imported routines:
  ReadMetadata

USE TableLib, ONLY : &
  !Imported definitions:
  Table, &
  TableCollection, &
  !Imported routines:
  TableNew, &
  TableGetNrows, &
  TableGetValue, &
  TableSetId, &
  TableSetTitle, &
  TableSetRowCol, &
  TableSetColHeader, &
  TableSetColUnit, &
  TableFillRow, &
  TableExport

USE Utilities, ONLY : &
  !Imported routines:
  GetUnit

USE IniLib, ONLY: &
  !Imported derived types:
  IniList, &
  !Imported routines:
  IniOpen, &
  IniClose, &
  IniReadReal, &
  IniReadString, &
  SectionIsPresent, &
  KeyIsPresent, &
  IniReadInt, &
  SubSectionIsPresent

USE Loglib, ONLY : &
  !Imported routines:
  Catch

USE StringManipulation, ONLY : &
  !Imported routines:
  ToString, &
  StringToFloat, &
  StringCompact, &
  StringToUpper, &
  StringTokenize, &
  StringToShort

USE GridLib, ONLY : &
  !Imported definitions:
  grid_real, &
  grid_integer

USE GeoLib, ONLY : &
  !Imported definitions:
  Coordinate, &
  !imported routines:
  DecodeEpsg

USE GridOperations, ONLY : &
  !Imported routines:
  GetIJ

USE Diversions, ONLY : &
    !Imported definition:
    Diversion

IMPLICIT NONE 

! Global (i.e. public) Declarations: 

TYPE Reservoir
	INTEGER(KIND = short)      :: id
	INTEGER(KIND = short)      :: rk !!runge-kutta order
	CHARACTER(LEN=100)         :: name
	CHARACTER(len=10)          :: typ !! on-stream off-stream by-pass
    TYPE (Coordinate)          :: xyz !!easting, northing and elevation in real world
	INTEGER(KIND = short)      :: r  !!cell row i
	INTEGER(KIND = short)      :: c  !!cell column j
    INTEGER(KIND = short)      :: rout  !!cell row  where off-stream pool outflow is discharged
	INTEGER(KIND = short)      :: cout  !!cell column where off-stream pool outflow is discharged
	REAL(KIND = float)         :: stage !!current stage [m]
	REAL(KIND = float)         :: stageMax !!maximum stage [m]
	REAL(KIND = float)         :: stageTarget !!follow a target (observed) stage [m]
	INTEGER(KIND = short)      :: typOut !!type ofoutflow: 1=free flow 2=target level
	INTEGER(KIND = short)      :: unit !!file unit  target stage
    INTEGER(KIND = short)      :: unitDischargeDownstream !!file unit of observed downstream discharge
    INTEGER(KIND = short)      :: unitDischargeDiverted !!file unit of observed diverted discharge
    LOGICAL                    :: dischargeDownstream  !!read observed downstream discharge when true
    LOGICAL                    :: dischargeDiverted  !!read observed diverted discharge when true
    INTEGER(KIND = short)      :: fileunit_out !!file unit for writing results
	TYPE(ObservationalNetwork) :: network ! pseudo station network of target file
    TYPE(ObservationalNetwork) :: networkDischargeDownstream ! pseudo station network of downstream discharge
    TYPE(ObservationalNetwork) :: networkDischargeDiverted ! pseudo station network of diverted discharge
	TYPE(DateTime)             :: tReadNewStage
    TYPE(DateTime)             :: tReadNewDischargeDownstream
    TYPE(DateTime)             :: tReadNewDischargeDiverted
	REAL(KIND = float)         :: Qout ! [m3/s]
    REAL(KIND = float)         :: Qout_off ! off-stream pool outflow [m3/s]
    REAL(KIND = float)         :: Pout_off ! off-stream pool outflow of previous time step [m3/s]
	REAL(KIND = float)         :: eFlow (365) !daily environmental flow [m3/s]
    REAL(KIND = float)         :: freeFlow ![m3/s] when Qin < freeFlow and stage <= freeFlowElevation Qout = Qin
    REAL(KIND = float)         :: freeFlowElevation !free flow elevation [m]
	TYPE(Table)                :: geometry !reservoir geometry and stage-discharge relationship
	TYPE(Table)                :: weir !used  by off-stream reservoirs and bypass channels
    INTEGER(KIND = short)      :: weirDOY (365) !weir function used for a given doy
    INTEGER(KIND = short)      :: geometryDOY (365) !geometry function used for a given doy
    REAL(KIND = float)         :: xout !! x coordinate where outflow from reservoir is discharged
    REAL(KIND = float)         :: yout !! y coordinate where outflow from reservoir is discharged
    LOGICAL                    :: highLevel !! true when reservoir is managed for high level
    REAL(KIND = float)         :: fullReservoirLevel !! full reservoir level (m)
    INTEGER                    :: qoutRule !!determines how to compute Qout 
    LOGICAL                    :: rising !!override Qout calculation only when Qin discharge is rising 
    TYPE (Diversion)           :: bypass !!diversion channel
    LOGICAL                    :: bypassIsPresent
    TYPE (Reservoir), POINTER  :: next  !dynamic list
    
END TYPE Reservoir

!Public routines
PUBLIC :: InitReservoirs
PUBLIC :: GetReservoirById
PUBLIC :: CountReservoirs
PUBLIC :: OutReservoirsInit
PUBLIC :: ReservoirSaveStatus

!Public variables
INTEGER (KIND = short) :: dtReservoir
TYPE (Reservoir), POINTER :: pools

!Local declarations:

!Local routines
PRIVATE :: ReservoirReadStatus
PRIVATE :: SetDailyArray
PRIVATE :: ReadGeometry
PRIVATE :: ReadWeir

!Local parameters
INTEGER, PARAMETER, PRIVATE  :: QIN = 1

!Local variables:
INTEGER (KIND = short), PRIVATE :: nReservoirs !total number of reservoirs
INTEGER (KIND = short), PRIVATE :: nReservoirsWithDiversion ! number of reservoirs with diversion


!=======
CONTAINS
!=======
! Define procedures contained in this module. 

  
!==============================================================================
!| Description:
!   Initialize reservoirs
SUBROUTINE InitReservoirs  &
!
(filename_ini,tbegin, mask, list)

IMPLICIT NONE

!arguments with intent in:
CHARACTER ( LEN = * ), INTENT(IN) :: filename_ini
TYPE (DateTime), INTENT(IN)       :: tbegin
TYPE (grid_integer), INTENT(IN)   :: mask


!arguments with intent out:
TYPE(Reservoir),POINTER,INTENT(out)   :: list !list of reservoirs

!local declarations
CHARACTER(LEN=1000)            :: filename
LOGICAL                        :: hotstart
TYPE(reservoir),POINTER        :: currentReservoir !points to current reservoir
INTEGER(KIND = short)          :: id_res !reservoir id
INTEGER(KIND = short)          :: k, i, j
TYPE(IniList)                  :: iniDB
INTEGER(KIND = short)          :: fileunit_geometry
INTEGER(KIND = short)          :: fileunit_weir
CHARACTER (LEN = 300)          :: string
LOGICAL                        :: error
INTEGER(KIND = short)          :: err_io
INTEGER(KIND = short), ALLOCATABLE :: doy (:)
INTEGER (KIND = short)         :: nDOY
INTEGER (KIND = short)         :: shortInt 
LOGICAL                       :: isString
!-------------------------------end of declaration-----------------------------

CALL Catch ('info', 'Reservoirs', 'Initializing reservoirs ')

!--------------------------------------------
!  open and read configuration file
!--------------------------------------------

CALL IniOpen (filename_ini, iniDB)

!--------------------------------------------
!  hot / cold start
!--------------------------------------------
IF ( KeyIsPresent ('path-hotstart', iniDB ) ) THEN
    hotstart = .TRUE.
ELSE
    hotstart = .FALSE.
END IF

!--------------------------------------------
!  allocate variables
!--------------------------------------------

nReservoirs =  IniReadInt ('nreservoirs', iniDB)

nReservoirsWithDiversion = 0 

!prepare list of reservoirs
 NULLIFY (list)
 DO k = 1, nReservoirs
   IF (.NOT. ASSOCIATED (list) ) THEN
     ALLOCATE (list)
     currentReservoir => list
   ELSE
     ALLOCATE (currentReservoir % next)
     currentReservoir => currentReservoir % next
   END IF
   
   !id
   currentReservoir % id = IniReadInt ('id', iniDB, section = ToString(k))
   
   !reservoir type: on-stream, off-stream
   currentReservoir % typ = &
       IniReadString ('type', iniDB, section = ToString(k))
   
   !name
   currentReservoir % name = &
       IniReadString ('name', iniDB, section = ToString(k))
   
   !coordinate
   currentReservoir % xyz % easting = &
       IniReadReal ('easting', iniDB, section = ToString(k))
   currentReservoir % xyz % northing = &
       IniReadReal ('northing', iniDB, section = ToString(k))
   currentReservoir % xyz % system = DecodeEpsg (IniReadInt ('epsg', iniDB))
   
   !local coordinate
   CALL GetIJ ( X = currentReservoir % xyz % easting, &
                Y = currentReservoir % xyz % northing, &
                grid = mask, i = currentReservoir % r, &
                j = currentReservoir % c )
   
   !runge-kutta order
   currentReservoir % rk = IniReadInt ('rk', iniDB, section = ToString(k))
   
    !initial stage
    IF (.NOT. hotstart) THEN
        currentReservoir % stage = &
          IniReadReal ('stage', iniDB, section = ToString(k))
    END IF
   
   SELECT CASE ( currentReservoir % typ )
   
   CASE ( 'on' ) !on-stream detention basin
       
       !target stage
       currentReservoir % eFlow = 0.
       IF (KeyIsPresent ('stage-target-file', iniDB, section = ToString(k) ) ) THEN
          string = IniReadString ('stage-target-file', iniDB, section = ToString(k))
          currentReservoir % unit = GetUnit ()
          OPEN ( unit=currentReservoir % unit , file=string)
          currentReservoir % typOut =  2
          currentReservoir % tReadNewStage = tbegin
          CALL ReadMetadata (currentReservoir % unit, currentReservoir % network )
          
          !environmental flow
          string = IniReadString ('e-flow', iniDB, section = ToString(k))
          
          currentReservoir % eFlow = SetDailyArray (string)
         
       ELSE
          currentReservoir % typOut =  1
       END IF
       
       !discharge downstream file
       IF (KeyIsPresent ('discharge-downstream-file', iniDB, &
                         section = ToString(k) ) ) THEN
          string = IniReadString ('discharge-downstream-file', iniDB, &
                                  section = ToString(k))
          currentReservoir % unitDischargeDownstream = GetUnit ()
          OPEN ( unit = currentReservoir % unitDischargeDownstream , &
                file = string)
          currentReservoir % tReadNewDischargeDownstream = tbegin
          CALL ReadMetadata (currentReservoir % unitDischargeDownstream, &
                             currentReservoir % networkDischargeDownstream )
          currentReservoir % dischargeDownstream = .TRUE.
       ELSE
          currentReservoir % dischargeDownstream = .FALSE.
       END IF
       
       
        !discharge diverted file
       IF (KeyIsPresent ('discharge-diverted-file', iniDB, &
                         section = ToString(k) ) ) THEN
          string = IniReadString ('discharge-diverted-file', iniDB, &
                                  section = ToString(k))
          currentReservoir % unitDischargeDiverted = GetUnit ()
          OPEN ( unit = currentReservoir % unitDischargeDiverted , &
                file = string)
          currentReservoir % tReadNewDischargeDiverted = tbegin
          CALL ReadMetadata (currentReservoir % unitDischargeDiverted, &
                             currentReservoir % networkDischargeDiverted )
          currentReservoir % dischargeDiverted = .TRUE.
       ELSE
          currentReservoir % dischargeDiverted = .FALSE.
       END IF
       
       
       !free flow
       IF  ( KeyIsPresent (key = 'free-flow', iniDB =  iniDB, &
                                       section = ToString(k) ) ) THEN 
            currentReservoir % freeFlow = IniReadReal ('free-flow', iniDB, &
                                                       section = ToString(k) )
       ELSE
            currentReservoir % freeFlow =  0.
       END IF
       
       !free flow elevation
       IF  ( KeyIsPresent (key = 'free-flow-elevation', iniDB =  iniDB, &
                                       section = ToString(k) ) ) THEN 
            currentReservoir % freeFlowElevation = &
                   IniReadReal ('free-flow-elevation', iniDB, &
                                        section = ToString(k) )
       ELSE
            currentReservoir % freeFlowElevation =  0.
       END IF
       
       !geometry
       CALL ReadGeometry (iniDB, k, currentReservoir )
       
       !manage high level options
       IF ( SubSectionIsPresent ('manage-high-level',ToString(k), iniDB ) ) THEN
            currentReservoir % highLevel = .TRUE.
           !read full reservoir level
           currentReservoir % fullReservoirLevel = &
                IniReadReal ( 'full-reservoir-level', iniDB, ToString(k), &
                              'manage-high-level')
           !read how to compute qout
           string = IniReadString ('qout', iniDB, ToString(k), &
                              'manage-high-level')
           SELECT CASE ( StringToUpper(string) )
           CASE ('QIN')
               currentReservoir % qoutRule = QIN
           END SELECT
           
           currentReservoir % rising =  &
               IniReadInt ( 'rising', iniDB, ToString(k), 'manage-high-level')
       ELSE
            currentReservoir % highLevel = .FALSE.
       END IF
       
       
   CASE ( 'off' ) !off-stream detention basin
       
       !maximum stage
       currentReservoir % stageMax = StringToFloat (string, error)
       currentReservoir % typOut =  1
       
       !geometry
       string = IniReadString ('geometry', iniDB, section = ToString(k))
       CALL TableNew (string, currentReservoir % geometry)
       
       !weir
       string = IniReadString ('weir', iniDB, section = ToString(k))
       CALL TableNew (string, currentReservoir % weir)
       
       !read x and y coordinate where outflow from reservoir is discharged
       currentReservoir % xout = &
           IniReadReal ('xout', iniDB, section = ToString(k))
       currentReservoir % yout = &
           IniReadReal ('yout', iniDB, section = ToString(k))
      
       CALL GetIJ ( X = currentReservoir % xout, &
                Y = currentReservoir % yout, &
                grid = mask, i = currentReservoir % rout, &
                j = currentReservoir % cout )
      
   CASE DEFAULT
       CALL Catch ('error', 'Reservoirs', &
            'unknown reservoir type: ', argument = currentReservoir % typ )
   END SELECT
   
   !diversion channel (optional)
   IF ( SubSectionIsPresent ('diversion',ToString(k), iniDB ) ) THEN
       
       currentReservoir % bypassIsPresent = .TRUE.
       
       nReservoirsWithDiversion = nReservoirsWithDiversion + 1
       
       !read x and y coordinate where outflow from reservoir is discharged
       currentReservoir % bypass % xout = &
           IniReadReal ('xout', iniDB, section = ToString(k), &
           subsection = 'diversion')
       currentReservoir % bypass % yout = &
           IniReadReal ('yout', iniDB, section = ToString(k), &
           subsection = 'diversion')
       
       !local coordinate in raster reference system
       CALL GetIJ ( X = currentReservoir % bypass % xout, &
                Y = currentReservoir % bypass % yout, &
                grid = mask, i = currentReservoir % bypass % rout, &
                j = currentReservoir % bypass % cout )
       
       !read weir data   
       CALL ReadWeir (iniDB, k, currentReservoir % bypass )
       
       !channel lenght
       currentReservoir % bypass % channelLenght = &
           IniReadReal ('channel-lenght', iniDB, section = ToString(k), &
                         subsection = 'diversion' )
       !channel slope
       currentReservoir % bypass % channelSlope = &
           IniReadReal ('channel-slope', iniDB, section = ToString(k), &
                        subsection = 'diversion' )
       !channel roughness coefficient
       currentReservoir % bypass % channelManning = &
           IniReadReal ('channel-manning', iniDB, section = ToString(k), &
                        subsection = 'diversion' )
       !channel section bottom width
       currentReservoir % bypass % channelWidth = &
           IniReadReal ('section-bottom-width', iniDB, section = ToString(k), &
                        subsection = 'diversion' )
       !channel section bank slope
       currentReservoir % bypass % channelBankSlope = &
           IniReadReal ('section-bank-slope', iniDB, section = ToString(k), &
                        subsection = 'diversion' )
       
       !environmental flow
       IF ( KeyIsPresent ('e-flow', iniDB, section = ToString(k),  &
              subsection = 'diversion' ) ) THEN
          string = IniReadString ('e-flow', iniDB, section = ToString(k),  &
            subsection = 'diversion' )
          currentReservoir % bypass % eFlow = SetDailyArray (string)
       ELSE !e-flow = 0.
         currentReservoir % bypass % eFlow = 0.
       END IF
   ELSE
       currentReservoir % bypassIsPresent = .FALSE.
   END IF
   
   NULLIFY (currentReservoir % next)
 END DO
 
 IF (hotstart) THEN
   string = IniReadString ('path-hotstart', iniDB)
   CALL ReservoirReadStatus (string)
 END IF
 
!--------------------------------------------
!  close configuration file
!--------------------------------------------

CALL IniClose (iniDb) 
 
RETURN
END SUBROUTINE InitReservoirs



!==============================================================================
!| Description:
!   initialise files for output
SUBROUTINE OutReservoirsInit &
  !
  (list, path_out)

IMPLICIT NONE


!arguments with intent(in):
TYPE(Reservoir), POINTER, INTENT(IN)   :: list !list of reservoirs
CHARACTER ( LEN = *), INTENT(IN) :: path_out

!local declarations:
TYPE(Reservoir), POINTER   :: currentReservoir !current reservoir in alist
CHARACTER (len = 100) :: string

!------------------------------end of declarations-----------------------------

currentReservoir => list
DO WHILE (ASSOCIATED(currentReservoir)) !loop trough all reservoirs
    
    string = ToString (currentReservoir % id)
    currentReservoir % fileunit_out = GetUnit ()
    OPEN(currentReservoir % fileunit_out,&
         file = path_out (1:LEN_TRIM(path_out))//'reservoir_'//&
               TRIM(string)//'.out')
     WRITE(currentReservoir % fileunit_out,'(a)') 'FEST: reservoir routing'
     WRITE(currentReservoir % fileunit_out,'(a,a)') &
         'reservoir name: ', currentReservoir % name &
        (1:LEN_TRIM(currentReservoir % name))
     IF ( currentReservoir % typ == 'off' ) THEN
        WRITE(currentReservoir % fileunit_out,'(a)')  'reservoir type: offstream'
     ELSE
        WRITE(currentReservoir % fileunit_out,'(a)')  'reservoir type: onstream'
     END IF
     IF ( currentReservoir % bypassIsPresent ) THEN
        WRITE(currentReservoir % fileunit_out,'(a)')  'with diversion: yes'
     ELSE
        WRITE(currentReservoir % fileunit_out,'(a)')  'with diversion: no'
     END IF
     WRITE(currentReservoir % fileunit_out,'(a,i4)') 'reservoir id: ', currentReservoir % id
     WRITE(currentReservoir % fileunit_out,*)
     WRITE(currentReservoir % fileunit_out,'(a)')'data'
     
     SELECT CASE ( currentReservoir % typ )
     CASE ( 'off' )
          IF ( currentReservoir % bypassIsPresent ) THEN
             WRITE(currentReservoir % fileunit_out,'(A)') 'DateTime h[m] Volume[m3] &
                Qupstream[m3/s] Qdownstream[m3/s] &
                 QinChannel[m3/s] QoutChannel[m3/s]'
          ELSE
             WRITE(currentReservoir % fileunit_out,'(A)') 'DateTime h[m] Volume[m3] &
                Qupstream[m3/s] Qdownstream[m3/s]'
          END IF
        
     CASE ( 'on' )
             WRITE(currentReservoir % fileunit_out,'(A)') 'DateTime h[m] Volume[m3] &
                 Qupstream[m3/s] Qdownstream[m3/s] &
                 QinChannel[m3/s] QoutChannel[m3/s]'
     END SELECT
    
     currentReservoir => currentReservoir % next
    
END DO

RETURN
END SUBROUTINE OutReservoirsInit
  
  
  
!==============================================================================
!| Description:
!   write results on files
SUBROUTINE OutReservoirs &
  !
  (list, time, Qin, Qout)

IMPLICIT NONE


!arguments with intent(in):
TYPE (Reservoir), POINTER, INTENT(IN) :: list !list of reservoirs
TYPE (DateTime),           INTENT(IN) :: time
TYPE (grid_real),          INTENT(IN) :: Qin
TYPE (grid_real),          INTENT(IN) :: Qout

!local declarations:
TYPE(Reservoir), POINTER   :: currentReservoir !current reservoir in alist
REAL (KIND = float)        :: volume !water stored in reservoir (m3)
REAL (KIND = float)        :: wse !water surface elevation 

!------------------------------end of declarations-----------------------------
timeString = time
currentReservoir => list
DO WHILE (ASSOCIATED(currentReservoir)) !loop trough all reservoirs
    wse = currentReservoir % stage
    CALL TableGetValue ( valueIn =  wse, tab = currentReservoir % geometry, &
                                     keyIn = 'h', keyOut ='volume', &
                                     match = 'linear', valueOut = volume, &
                                     bound = 'extendlinear' )
    SELECT CASE ( currentReservoir % typ)
    CASE ( 'off' )
        IF ( currentReservoir % bypassIsPresent ) THEN
            WRITE(currentReservoir % fileunit_out,'(a,6(" ",e14.7) )') timeString,&
              currentReservoir % stage, volume, &
              Qin % mat(currentReservoir % r, currentReservoir % c),&
              Qout % mat(currentReservoir % r, currentReservoir % c), &
              currentReservoir % bypass % QinChannel, &
              currentReservoir % bypass % QoutChannel
        ELSE
             WRITE(currentReservoir % fileunit_out,'(a,4(" ",e14.7) ) ') timeString,&
               currentReservoir % stage, volume, &
               Qin % mat(currentReservoir % r, currentReservoir % c),&
               Qout % mat(currentReservoir % r, currentReservoir % c)
        END IF 
    CASE ( 'on' )
        IF ( currentReservoir % bypassIsPresent ) THEN
            WRITE(currentReservoir % fileunit_out,'(a,6(" ",e14.7) )') timeString,&
              currentReservoir % stage, volume, &
              Qin % mat(currentReservoir % r, currentReservoir % c),&
              Qout % mat(currentReservoir % r, currentReservoir % c), &
              currentReservoir % bypass % QinChannel, &
              currentReservoir % bypass % QoutChannel
        ELSE
             WRITE(currentReservoir % fileunit_out,'(a,6(" ",e14.7) )') timeString,&
               currentReservoir % stage, volume, &
               Qin % mat(currentReservoir % r, currentReservoir % c),&
               Qout % mat(currentReservoir % r, currentReservoir % c), &
               0.0, 0.0
        END IF
    END SELECT
    
     currentReservoir => currentReservoir % next
    
END DO

RETURN
END SUBROUTINE OutReservoirs 



!==============================================================================
!| Description:
!   Save reservoirs status into file
SUBROUTINE ReservoirSaveStatus &
  !
  ( pathOut, time )

IMPLICIT NONE


!arguments with intent(in):
CHARACTER ( LEN = *), INTENT (IN) :: pathOut
TYPE (DateTime), OPTIONAL, INTENT (IN) :: time

! local declarations:
TYPE (Table) :: tab
CHARACTER (LEN = 300) :: string
CHARACTER (LEN = 100), ALLOCATABLE :: row (:)
INTEGER (KIND = short) :: i
TYPE (Reservoir), POINTER :: currentReservoir !points to current reservoir
!-------------------------------end of declaration-----------------------------

!create new table for reservoir stage
CALL TableNew ( tab )

!populate table
string = 'reservoir stage'
CALL TableSetId ( tab, string)

IF ( PRESENT (time) ) THEN
  timeString = time
  string = 'reservoir stage at: ' // timeString
ELSE
  string = 'reservoir stage at the end of simulation' 
END IF
CALL TableSetTitle ( tab, string)

!Allocate variables
CALL TableSetRowCol ( tab, nReservoirs, 2 )
IF ( ALLOCATED ( row ) ) THEN
    DEALLOCATE ( row )
END IF
ALLOCATE ( row (2) )

!set column header and unit
CALL TableSetColHeader (tab, 1, 'id')
CALL TableSetColHeader (tab, 2, 'stage')

CALL TableSetColUnit (tab, 1, '-')
CALL TableSetColUnit (tab, 2, 'm')

!fill in rows
currentReservoir => pools

DO i = 1, nReservoirs
     !id
     row (1) = ToString (currentReservoir % id)
     !stage
     row (2) = ToString (currentReservoir % stage)
     
     currentReservoir => currentReservoir % next
     
     CALL TableFillRow (tab, i, row)
END DO

IF (PRESENT(time)) THEN
	timeString = time
	timeString = timeString (1:19) // '_'
	timeString (14:14) = '-'
	timeString (17:17) = '-'
		
ELSE
	timeString = '                    '
END IF

string = TRIM(pathOut) // TRIM(timeString) // 'reservoirs.tab'

CALL TableExport (tab, string )


!create new table for diverted discharge
CALL TableNew ( tab )

!populate table
string = 'diverted discharge'
CALL TableSetId ( tab, string)

IF ( PRESENT (time) ) THEN
  timeString = time
  string = 'diverted discharge at: ' // timeString
ELSE
  string = 'diverted discharge at the end of simulation' 
END IF
CALL TableSetTitle ( tab, string)

!Allocate variables
CALL TableSetRowCol ( tab, nReservoirsWithDiversion, 3 )
IF ( ALLOCATED ( row ) ) THEN
    DEALLOCATE ( row )
END IF
ALLOCATE ( row (3) )

!set column header and unit
CALL TableSetColHeader (tab, 1, 'id')
CALL TableSetColHeader (tab, 2, 'Qin')
CALL TableSetColHeader (tab, 3, 'Qout')

CALL TableSetColUnit (tab, 1, '-')
CALL TableSetColUnit (tab, 2, 'm3/s')
CALL TableSetColUnit (tab, 3, 'm3/s')

!fill in rows
currentReservoir => pools

DO i = 1, nReservoirs
    IF ( currentReservoir % bypassIsPresent ) THEN
     !id
     row (1) = ToString (currentReservoir % id)
     !Qin
     row (2) = ToString (currentReservoir % bypass % QinChannel)
     !Qout
     row (3) = ToString (currentReservoir % bypass % QoutChannel)
     
     currentReservoir => currentReservoir % next
     
     CALL TableFillRow (tab, i, row)
    END IF
END DO

IF (PRESENT(time)) THEN
	timeString = time
	timeString = timeString (1:19) // '_'
	timeString (14:14) = '-'
	timeString (17:17) = '-'
		
ELSE
	timeString = '                    '
END IF

string = TRIM(pathOut) // TRIM(timeString) // 'reservoirs.tab'

CALL TableExport (tab, string, append = .TRUE. )

RETURN
END SUBROUTINE ReservoirSaveStatus
  
  
  
!==============================================================================
!| Description:
!   Read reservoir stage and diversion discharge from file
SUBROUTINE ReservoirReadStatus &
  !
  (filename)

IMPLICIT NONE


!arguments with intent(in):
CHARACTER ( LEN = *), INTENT (IN) :: filename

! local declarations:
INTEGER (KIND = short) :: i
CHARACTER (LEN = 300)  :: string
TYPE (TableCollection) :: tabs
TYPE (Reservoir), POINTER :: currentReservoir !points to current reservoir

!------------------------------end of declarations-----------------------------

!read tables
CALL TableNew ( filename, tabs )

!set reservoir stage
currentReservoir => pools
DO i = 1, nReservoirs
    string = TRIM ( ToString ( currentReservoir % id ) )
    
    CALL TableGetValue (  valueIn = string, &
                          tables = tabs, &
                          id = 'reservoir stage', &
                          keyIn = 'id', &
                          keyOut = 'stage', &
                          valueOut = currentReservoir % stage )

    currentReservoir => currentReservoir % next
END DO 

!set diversion discharge
IF ( nReservoirsWithDiversion > 0 ) THEN
    currentReservoir => pools
    DO i = 1, nReservoirs
        IF ( currentReservoir % bypassIsPresent ) THEN
            string = TRIM ( ToString ( currentReservoir % id ) )
    
            CALL TableGetValue (  valueIn = string, &
                        tables = tabs, &
                        id = 'diverted discharge', &
                        keyIn = 'id', &
                        keyOut = 'Qin', &
                        valueOut = currentReservoir % bypass % QinChannel )
            
             CALL TableGetValue (  valueIn = string, &
                        tables = tabs, &
                        id = 'diverted discharge', &
                        keyIn = 'id', &
                        keyOut = 'Qout', &
                        valueOut = currentReservoir % bypass % QoutChannel )
        END IF
        currentReservoir => currentReservoir % next
    END DO 
END IF

RETURN
END SUBROUTINE ReservoirReadStatus  
  
!==============================================================================
!| Description:
!   given a list of reservoirs, returns a pointer to one reservoir by id
FUNCTION GetReservoirById &
( list, id) &
RESULT (p)

IMPLICIT NONE 

!Arguments with intent in:
TYPE(Reservoir), POINTER, INTENT(IN) :: list !list of reservoirs
INTEGER(KIND = short), INTENT(IN)    :: id

!Arguments with intent out:
TYPE(Reservoir), POINTER :: p !pointer to reservoir

!local arguments:
LOGICAL :: found

!------------end of declaration------------------------------------------------

!loop througout list of reservoirs searching for id
p => list
found = .false.
DO WHILE (ASSOCIATED(p)) 
  IF (p % id == id) THEN
    found = .TRUE.
    EXIT
  ELSE
    p => p % next
  END IF
ENDDO

IF (.NOT. found ) THEN
  CALL Catch ('error', 'Reservoirs', 'reservoir not found by &
                       function GetReservoirById: ', argument = ToString(id))
END IF

RETURN
END FUNCTION GetReservoirById
  


!==============================================================================
!| Description:
!   count number of reservoirs in a list
FUNCTION CountReservoirs &
( list) &
RESULT (c)

IMPLICIT NONE 

!Arguments with intent in:
TYPE(Reservoir), POINTER, INTENT(IN) :: list !list of reservoirs

!Arguments with intent out:
INTEGER(KIND = short) :: c

!local arguments:
TYPE(Reservoir), POINTER :: p !pointer to reservoir

!------------end of declaration------------------------------------------------

!loop througout list of reservoirs searching for id
p => list
c = 0
DO WHILE (ASSOCIATED(p)) 
  c = c + 1
  p => p % next
ENDDO

RETURN
END FUNCTION CountReservoirs


!==============================================================================
!| Description:
!   populate  array of daily values
FUNCTION SetDailyArray &
!
( value ) &
!
RESULT ( array )

IMPLICIT NONE

!Arguments with intent(in):
CHARACTER (LEN = *) :: value


!Local declarations:
REAL (KIND = float) :: array (365)
REAL (KIND = float) :: scalar
LOGICAL :: error
TYPE (Table) :: valueTable
INTEGER :: i
REAL (KIND = float) :: doyStart, doyStop

!----------------------------end of declarations-------------------------------

!check that value is a number
scalar = StringToFloat (value, error)
IF  ( error ) THEN !value changes in time
    CALL TableNew (value, valueTable)
    array = 0.
    
    DO i = 1, TableGetNrows (valueTable)
        CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='doy-start', &
                             match = 'exact', &
                             valueOut = doyStart )
      
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='doy-stop', &
                             match = 'exact', &
                             valueOut = doyStop )
       
       CALL TableGetValue ( valueIn = REAL (i), &
                             tab = valueTable, &
                             keyIn = 'count', &
                             keyOut ='value', &
                             match = 'exact', &
                             valueOut = scalar )

       array ( INT (doyStart) : INT (doyStop) ) = scalar
        
    END DO
    
ELSE !no error, value is a scalar
    array = scalar
END IF

RETURN
END FUNCTION SetDailyArray


!==============================================================================
!| Description:
!   read geometry table
SUBROUTINE ReadGeometry &
  !
  (iniDB, k, reserv)

IMPLICIT NONE

!Arguments with intent(in):
TYPE(IniList), INTENT (IN) :: iniDB
INTEGER (KIND = short), INTENT (IN) :: k

!Arguemnts with intent (inout):
TYPE(Reservoir), POINTER, INTENT(INOUT) :: reserv 


!local declarations
CHARACTER (LEN = 300)          :: string
INTEGER (KIND = short)         :: nDOY
INTEGER (KIND = short)         :: shortInt 
INTEGER (KIND = short)         :: i, j
INTEGER(KIND = short), ALLOCATABLE :: doy (:)
LOGICAL                       :: isString

!---------------------------end of declarations--------------------------------

string = IniReadString ('geometry', iniDB, section = ToString(k))
CALL TableNew (string, reserv % geometry)

nDOY = reserv % geometry % noCols - 3
ALLOCATE ( doy ( nDOY ) )

j = 0
DO i = 1, reserv % geometry % noCols
    shortInt = StringToShort ( reserv % geometry % col ( i ) % header, isString )
    
    IF ( .NOT. isString ) THEN !number detected
        j = j + 1
        doy ( j ) = shortInt
    END IF
    
END DO

!set  when (doy day of year) geometry table changes
reserv % geometryDOY (:) = MAXVAL (doy)
       !
DO i = 1, nDOY
    DO j = doy (i), 365
        reserv % geometryDOY (j) = doy (i)
    END DO
END DO
       
DEALLOCATE (doy)

RETURN
END SUBROUTINE ReadGeometry

!==============================================================================
!| Description:
!   read weir table
SUBROUTINE ReadWeir &
  !
  (iniDB, k, div)

IMPLICIT NONE

!Arguments with intent(in):
TYPE(IniList), INTENT (IN) :: iniDB
INTEGER (KIND = short), INTENT (IN) :: k

!Arguemnts with intent (inout):
TYPE(Diversion), INTENT(INOUT) :: div 


!local declarations
CHARACTER (LEN = 300)          :: string
INTEGER (KIND = short)         :: nDOY
INTEGER (KIND = short)         :: shortInt 
INTEGER (KIND = short)         :: i, j
INTEGER(KIND = short), ALLOCATABLE :: doy (:)
LOGICAL                       :: isString

!---------------------------end of declarations--------------------------------

string = IniReadString ('weir', iniDB, section = ToString(k), &
                         subsection = 'diversion')
CALL TableNew (string, div % weir)

nDOY = div % weir % noCols - 1
ALLOCATE ( doy ( nDOY ) )

j = 0
DO i = 1, div % weir % noCols
    shortInt = StringToShort ( div % weir % col ( i ) % header, isString )
    
    IF ( .NOT. isString ) THEN !number detected
        j = j + 1
        doy ( j ) = shortInt
    END IF
    
END DO

!set  when (doy day of year) geometry table changes
div % weirDOY (:) = MAXVAL (doy)
       !
DO i = 1, nDOY
    DO j = doy (i), 365
        div % weirDOY (j) = doy (i)
    END DO
END DO
       
DEALLOCATE (doy)

RETURN
END SUBROUTINE ReadWeir

END MODULE Reservoirs